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f»0 , We derive a simple analytical expression for the two-body force in a sub-class of MOND-like 

theories and make testable predictions in the modification to the two-body orbital period, shape, 
and precession rate, and escape speed etc. We demonstrate the applications of the modified Kepler's 
law in the timing of satellite orbits around the Milky Way, and checking the feasibility of MOND 
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in the orbit of Large Magellanic Cloud, the M31 galaxy, and the merging Bullet Clusters. 
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MOND 
appears to be consistent with satellite orbits although with a tight margin. Our results on two-bodies 
are also generalized to restricted three-body, many-body problems, rings and shells. 

PACS numbers: 98.10.+Z, 95.35.+d, 98.62.Dm, 95.30.Sf 
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The Kepler's law, or the full analytical solution to the two-body problem is perhaps the most powerful prediction 
O ■ of Newton's gravitational law, iFd = I — F2I = Gm j m2 for the forces between any two point masses mi, ni2 orbiting 

-^ , each other with a separation r\i. For bound orbits the Kepler's law predicts the relations between masses and the 
C0 . r-73- 



INTRODUCTION 



orbital period (T oc V 7^+^177 )• ^ ne total mass of any binary or merging cosmological bodies can be constrained 
by the fact that first a binary must be bound, and second the binary's period has to be smaller or comparable to the 
age of the universe. These simple conditions of boundness and timing have powerful applications in constraining the 
unknown physics of the widely-speculated cosmological dark matter. 

For example, if the Milky Way had nothing except its baryonic mass mi ~5x 10 10 solar masses, then the specific 

L> J gravitational potential ^r 1 would be too shallow to bound the fast moving satellites of the Milky Way, e.g., the LMC 
and Sgr with velocity about 300-380 km/s at 50 kpc and 16 kpc from our center well exceeds their local escape speed 
of about 100-150km/s. To explain the long tidal tails around these satellites, one must deepen the potential, e.g., by 
adding a dark matter halo, so that these satellites can be on bound orbits and have a few pericenter passages to create 
the tidal tails. There are constraints coming from being bound and constraint Another example is the colliding Bullet 

7-H . clusters, which shows enormous relative velocity in the shocked X-ray gas ~ 4000 km/s; in comparison if the merging 
clusters were modeled as two baryonic gas clouds of total mass < 10 14 M Q falling for the first time from infinity to the 



present separation of 400-700 kpc, one expects at most a Keplerian escape speed J — (™i+ m 2) <• I500km/s. It has 
;_J ■ been suggested that even the gravity of cosmologically reasonable amount of dark matter is barely enough, and an 
attractive fifth force between dark matter might be necessary to explain the fast speed of the bullet in a Newtonian 
gravity [7, 10, 21]. 

Likewise the Kepler's law allows us to time the motions of celestial objects. For example, in the limit that the 
local group can be approximated as two points masses: the Milky Way and M31 galaxy, Kepler's law can be used to 
estimate their total mass mi + m^ ~ (2 — 3) x 1O 12 M , one order of magnitude larger than the baryonic mass, hence 
one can argue the existence of dark halos in these two galaxies under Newton's gravity if 7*12 is comparable to the 
semi-axis of the orbit, and the half-period T/2 ~ 14 Gyrs so that the two has enough time to move apart and turn 
around in a Hubble time. 

Nevertheless there is another way to make galaxy potential deeper without introducing unknown matter. The 
Modified Newtonian Dynamics (MOND) predicts an enhancement of the gravitational coupling constant when the 

gravity drops below o,q ~ 1.2 x 10~ 10 m/s 2 , G ~ G 1, -nJV . This produces the effect of a dark halo without 

L '*U max 

actually invoking real dark matter [3]. It is interesting to derive the modified Kepler's law for the MOND theories, 
and contrast their predictions on satellite orbits with Newton's gravity. Unfortunately, the prediction of the MOND 



theories generally invokes solving a modified Poisson equation V • — —^ — ^ m,i5(r — r,) numerically over all 3- 
dimensional space for the gravity g before obtaining the forces on the two-bodies. The problems become even harder 
if there are N > 2 bodies. To study the orbit, one would need an N-body code to integrate the orbits step by step 
and solve for the modified Poisson equation each step. To date there have been no realistic simulations of the orbits 
of two point- masses in MOND. 

Here we show that in some versions of MOND, the virial is a fully analytical expression for arbitrary matter 
densities. In the limit that we can neglect the sizes of the bodies, one can invert the virial 7-12-F12 to find the forces 
-F12 between two bodies. This completely by-passed the Modified Poisson equation. The advantage of this analytical 
result, rigorous in limited versions of the MOND theory, is to facilitate the estimation of the modification effects 
of MOND-like theories. E.g., Previous calculations of the timing of the encounter of the Bullet Cluster in MOND 
are based purely on numerical simulations [1, 9, 17], which made the essential physics somewhat more obscure. Our 
result allows one to gain similar level of analytical intuition as enjoyed in standard Dark Halo theories in Newtonian 
gravity, albeit the analytics can never substitute more rigorous and realistic numerical cosmological calculations and 
simulations in all these theories. 

It is worth recalling two interesting features in Newton's gravity: the pericentcr of the bodies has no precession 
because there is no distinction of a radial oscillation period and a period to turn 360 degrees in angular direction. 
The virial of the two-body system W = X)i=i r i ' ^i equals to — Gmini2 ; which applies at any time, whether the orbit 
is elliptical or hyperbolic. We will discuss how these properties are modified in MOND. Note that the concept of 
virial applies in all classical gravity theories, and is more general than the Newtonian virial theorem, and one can 
speak of the instantaneous virial without any assumption of time-averaging or equilibrium state of the objects. 

The outline of the present work is as follows: in §2 we give the actions for two specific versions theories for MOND, 
and show the simple expression for the virial in these theories. In §3 derives the two-body force from the virial, 
and give the equation of motion. §4 generalizes the results to many-body, rings and shells. In §5 we illustrate the 
applications of the analytical results in calculating the orbits of the Local Group objects. We summarize in §6. 

HOW TO CALCULATE THE VIRIAL IN NEWTONIAN AND MONDIAN GRAVITY 

For calculability, we adopt the multi- field version of MOND according to Bekenstein [2, 3] or the recent Qusai-linear 
MOND (QMOND) version of Milgrom[15]. In the low-speed weak field limit, particles move under their Newtonian 
potential &n and the MONDian scalar field potential $ s (which plays the role of Dark Matter potential). These two 
potentials are given by following Poisson equations for an N-body system 
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V 2 $ w = 4ttGp , V 2 $ s = V • (^ Q mcwdV$ w ), vqmond = I L ^ LL J , m QMOND (1) 

V 2 $ w - 4ttGp , V • (pmondV^s) = V 2 $at, ^mond = \—^) » in MOND, (2) 

where the uqmon D-iw^ction or the /iMCWD-function leads to a dcep-MOND effect, and the matter density, 

N N 

consists of N softened particles, each with the Lagrangian density 
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where rrii/V(r — rj(i), b) is a spherical density profile with a softening radius b of the particle center rj(t), e.g., a so- 
called Hernquist density profile. [26] The softened particles ensure that we do not over-generalize the above theories 
to situations of strong gravity. 

The above modified Poisson equations can be obtained self-consistently from minimizing the following action with 
respect to &n or $ s : 
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Smond = / dtdr 3 l^L 



N 

N 



V$/v| 2 . |V$ S 



8nG 12nGa 

The equation of motion for both versions of MOND is derived by variation of the action S with rj, which gives 
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dt dt adt 2 dvi 

Here the particles are coupled to the total potential 

$ = $at + $ s (8) 

which is consisted of two parts, the Newtonian part $tv and the scalar field part $ s . The gradient of the scalar field 
— V$ s in both versions of MOND plays the effective role of the dark matter gravity [23, 25], which is predicted here 
of the amplitude (Gao^2 i m,i) 1 / 2 /r at distances r far from the N particles, i.e., a test particle at a large distance r 
would accelerate with 

9($at + $ s ) ^ r ,v^ n± r /■n^ 

~^^ — ^«-L Gm v-(L Ga ° m ^v (9) 
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To see that our formulation indeed recovers the Newtonian and deep-MOND limits in a static universe, we note 
that at small radii from a mass m,, a test particle's acceleration g goes as Grriir~ 2 , and at large radii the acceleration 
g goes as (Gaoirii) ?r~ l . With a bit of algebra one can show the equivalent /i-function in MOND is given by p{g) = 



M. =v -i = i_ 
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[24]. 
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if one can assume the distribution of the matter and gravity has spherical symmetry 



APPLICATION TO TWO-BODY PROBLEM IN GENERAL 

Zhao & Famaey [25] found that, in the absence of cosmic expansion, the virial W and its Newtonian counterpart 
Wn satisfy a simple relation 
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\W\= dr 3 pr-(V$ N + V$ s ) = \W N \ + -JGao(y2m l Y, \W N \= dr 3 pr-V$ N (10) 

Jo •> Jo 

which applies to both QMOND and multi-field MOND for an isolated matter distribution in any geometry. 

We now apply this to a general two-body system. Following Milgrom[14], we argue that the total virial can be 
broken apart into the internal part and orbital part. In the limit the bodies are separated with distances much bigger 
than their sizes so that the mutual gravity is much weaker than the internal gravity in the vicinity of each (compact) 
mass, we can neglect the external field when applying Poisson's equation inside each body, each body can be treated 
as isolated system, hence each satisfies its own virial theorem. 

From this we can estimate the internal virial, 

\Wi\ = j dr 3 pv • V*am + l^Gmfao. (11) 

Subtracting off the internal virial J^ Wi from the total virial W , we find the interaction virial r\iF\i is given by 
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ri2 • P12 = ^^ + ^G(m 1+ m 2 )3 ao - ]T =-jGm 3 a . (12) 
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where we have used the fact that the total Newtonian potential energy subtracting the Newtonian potential energy 
of each body yields just the mutual potential energy. As far as the particles have very small sizes compared to their 
separation, we have ri2 • -F12 = ^yiFyi^ where ri2 has a negligible spread in distance. So we find that the mutual force 
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Clearly the force is Newtonian at short distance, and at large distance tends to a MONDian 1/r force with a non-trivial 
normalization. 

The accelerations of the bodies ^^ and — £12., or the gravity F 12 between bodies are actually independent of their 
relative velocities and history, hence the expression is general for the mutual (MONDian) gravity between two (widely- 
separated compact) bodies. The above expression for force is rigourous if the sizes of bodies are smaller than their 
distance. In reality the (spherical) particles can have a finite size, say b, so that the distance between two bodies have 
a distribution of width ~ b instead of a single value |ri2 1- Extended bodies also introduce new effects, such as studied 
in [6] in the context of the Birkhoff theorem. Nevertheless one could adapt the formulae to approximate the effect of of 
non-spherical bodies or softened particles. Such a formula is given in Appendix for two axisymmetrical particles with 
a bulge of length scale b and disk scale length k. Especially in calculating the force between two spherical particles, 
one replaces r\i — > rri + b in eq. 13. In all cases we ensure a construction of MOND two-body force that is conserving 
momentum, energy and angular momentum of the whole system. 

Note that the two-body force here is far from obvious. A naive application of MOND could often lead to incorrect 

answer, e.g., iq = Gm r \ m2 + m 2 J Gm r .l a ° - or F 2 = Gm r \ m2 + m 1 J G ^ a ° -, which would violate momentum-conservation 
[1] . Also our two-body result does not hold rigorous if adopting the Bekenstein-Milgrom (BM84) theories of MOND 
where the scalar field is not used, e.g., the calculations in [6]. In general the forces must be computed numerically by 
first solving the Poisson equations and then integrating the force over the volume of the (extended) body concerned. 
Our analytical result here helps to calibrate numerical grid or boundary effects in a numerical code. 

Two-body equation of motion in cosmological background 

We are almost ready to apply our derived force to the two-body problem except that we have to consider the 
effect of the Hubble expansion, which is non-negligible for any timing arguments. Considering the expansion of the 
background universe a(t) so that the relative distance of particles 1 and 2 in proper coordinates iq — r 2 = (xi — x 2 )a(£), 
we find the equation of motion is given by 

dt dt adt 2 \TTiiJ 

Note here the factional term % 4z± i n equations in co- moving coordinates has canceled itself when the equation is 
written for the proper coordinates r^. Approximate the scale factor as a(t) = {t/to) n , then the cosmological term 
is n(n — 1)£~ 2 , and is zero if n — (static) or n = 1 (empty open universe). For LCDM parameters, the Hubble 
parameter ^ = 14 ^ \/0-667 + 0.333a~ 3 , we find the following approximation, a(t) oc (e y — l) 2 / 3 , where y = i/llGyr 

and j^ = 3(11 q rs )i eV (^ eV — l)( e!y — I) -2 so that we expect a nearly constant £^ — 9(11 q TS yi at late times and in 

the future, and -^A- = -53 at earlier times when a(t) ex t 2 / 3 . Generally the cosmological term is an attractive force 
in matter-domination, and repulsive in vaccum domination. For the latter part, we shall consider only a universe 

dominated by a pure vaccum energy density p c — g^g, the cosmological term 

gC = ^ ri2 = T ri2 ' T= {-^-) (15) 

acts as a static repulsive potential force for an exponentially growing scale factor a(t) = exp(£-\/87rG / o c /3) for positive 
p c . We adopt a vacuum universe with the cosmic e- folding time r w 11 Gyrs. 
The two-body relative acceleration 
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MANY-BODY PROBLEM IN MOND 



Consider a 2D or 3D symmetric distribution of N > 1 identical particles of mass rrii — m/N plus a central particle 
of mass M. As long as the symmetry guarantees that the forces on the particles m, is pointing to the central mass 
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FIG. 1: shown is an example of A^ = 6 self-gravitating particles each with mass m* = M/(50N) initially on a ring of radius 16 
kpc moving 280 km/s tangentially around the Milky Way (red point at origin M = 5 x 10 10 solar masses); the arrows indicate 
the force towards the central mass; one can see the precession of the pericenter and the stretching of the ring for the past 1 
Gyrs. The right panel shows a semi-log plot of the Milky Way's potential 2V(r) (in units of square km/s) for a tracer particle 
(solid) and N = 6 massive particles on a ring of mass A-//2 (dashed line) as function of its distance r, where we have labeled 
the location of the cut off radius, where 2V(r cut ) ~ 0; This allows some stars to escape. In contrast stars are unable to escape 
from a logarithmically divergent potential of an isolated MOND galaxy in an empty universe. All length units are in kpc. 



M, and that particle M does not experience any force, we can use the same virial theorem to obtain a relation of the 
acceleration g of each particle at |r^| = r, and the Newtonian acceleration gjq and the cosmological acceleration gc as 
follows 
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(Gm, 3 «o) 1/2 



(18) 



Note that the first zero stands for the zero-virial acting on the central mass M. Also note that any non-symmetry or 
finite sizes would prevent us inverting the virial for the force because the force on each particle would not all be of the 
same amplitude and pointing radially. Such symmetric configurations are generally contrived with no counterparts in 
reality, except for perhaps N = 2 for binary systems, and N ~ oo for astronomical rings and shells, which are fairly 
common due to mergers of galaxies, which can form polar rings around spiral galaxies, and shells around elliptical 
galaxies. 

Making a correction for the expansion of the vacuum only universe of a constant cosmic density p c , and taking into 
account of the finite size b, we get 
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(19) 



where j is the specific angular momentum of the system, which can be solved from the pericenter radius r per i , where 
dr/dt = 0. Here |j is a dimensionless parameter depending on the mass ratios and the detailed geometry of the 
system and a at = win I'm is a geometrical factor taking into account of the addition of Newtonian forces among the 
N particles. 
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where 2^ = | 



1 



M' 1 



'N~ 



(M+m)3/2 



Note that the N = 1 case is special because the particle M would move around 

the center of the mass to balance with the single particle m conserve total momentum, while the particle M coincides 
with the center of mass for N > 1 symmetrically distributed particles of identical mass m/N . In the appendix we 
give two expressions which shows the asymptotic behavior more clearly in cases with extreme mass ratio. We find 
M -¥ M for small mass ratio m/M — > 0, and M — > (1 - A r_1 / 2 ) 2 (4/9)m for large mass ratio m/M — > oo. If N = 1, 
we have M —>■ M or m for extreme mass ratios. It can be shown that 0.6 (M + m) < M < (M + m) for any mass 
ratios and N. 



For N point masses distributed on regular polygon, we find 
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for N = 6 etc. etc.. By 



specifically a 2 = § for TV = 2, a 3 = ^ for TV = 3, a 4 = ^^^ for N = 4, a 6 

choosing a large TV we can essentially model a circular "ring" of self-gravitating particles. A similar calculation can 
be done for a^ if the points are distributed on regular polyhedra; the expressions are given Appendix C. For large 
N the configuration is a fair approximation to a spherical "shell" distribution of self-gravitating particles. 

The above analytical results reveal several interesting distinctions of the MOND force with Newtonian force, (i) 
On a ring or a polygon with A~ — > oo, the Grh/r 2 term in the Newtonian force diverges, and V Ga$M /r term in 
the MOND part approaches a common constant rotation curve because the MONDian M 1 / 2 /M 1 / 2 term depends 
on the mass ratio m/M and linearly on A" 1 ' 2 , while in Newtonian the rh/m term depends purely non-linearly on 
the geometry parameter N — 1. (ii) The MONDian force gives the same flat rotation curve while the Newtonian 
force produces Keplerian rotation curves of different amplitudes whether the N particles form a polygon/ring or a 
polyhedron/shell, (iii) The MONDian force does produce precession while the Newtonian force does not. 



Effective potential, orbital period and precession 



The equation of motion eq. 19 can be reduced to purely radial motion as 

fr = _d VM -3 2 r 2 G(M + m N ) 

dt 2 dr K h K >~ 2r 2 2r 2 r + b 



(GMa ) 1/2 Hr + b) 



(23) 



where we absorb the centrifugal force into an effective potential V(r), j is the specific angular momentum of the 
system. We can integrate the equation of motion using energy conservation to get 
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(24) 



For non-radial motion with an angular momentum barrier the effective potential relates also the pericenter with 
the apocenter via V(r peri ) — V(r apo ) = E, where dr/dt = 0. The time from the pericenter to the apocenter and then 
turn back to the pericenter is given by 



2dr 

\dr/dt\ ' 



(25) 



Each of the particles on a polygon/ring will make a rosette orbit, which can precess while keeping the configuration 
self-similar; the whole pattern resembles a closing/opening shutter. The backward precession angle per orbit is 
determined by [4] 
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In the case that the N particles are on a circular orbit of speed v C i r , we have a potential of N-fold rotational symmetry 
with a pattern rotation angular speed u> pa ttem — v C i r (r)/r and the rotation curve v c i r (r) is given by 



Vcir(r) = 



r 2 G(M + m N ) r 



(r + b) 



■(GMao) 1 ' 2 



-,1/2 



(27) 



where the angular momentum j is found by requiring the effective potential satisfies dV/dr = at radius r = r apo = 
Tperi- Clearly the rotation curve v c ir(r) is flat at large radii if there is no cosmological background. The pattern 
rotation period or pattern speed of the potential is especially useful for modeling a self-consistent m = N = 2 bar 
potential in MOND. Note that the precession here does not change the orbital plane fixed by the constant angular 



momentum vector, only the direction of the pericenters, in some sense similar to mercury's pericenter precession 
due to GR corrections to the Kepler force of two bodies. Any precession of the orbital plane would be the effect of 
non-spherical shapes of the bodies, which we do not go into in this paper. 

Fig. 1 gives an illustration of the many-body problem in MOND, where the orbital parameters are inspired by 
those of the Sgr satellite, whose mass is broken into a self-gravitating tidal debris on a polar plane around the Milky 
Way. We approximate the Sgr mass by N = 6 massive particles on a circle of 16 kpc around the central mass M, 
representing the Milky Way. One can clearly see the precession of the pericenter, the rotation/expansion/shrinking 
of the initial ring. 

For radial motion with j = r per i = 0, the speed that the particles cross the origin is given by yj2(y(r apo ) — V(0)). 
E.g., the N particles can move radially on a set of self-similiar regular polyhedron (an approximation to more realistic 
shells). For the radial motion, there is a maximum (or an edge) of the effective potential due to the repulsive force 
gc balancing the Newtonian and MONDian inward force, 
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Neglecting the r~ 2 term due to Newtonian force, and neglecting the finite size 6, we get the edge of the potential is at 

Tedge = (GMa ) 1/4 T. (29) 

The particles in the MONDian system will not be bound if the TV particles have a specific energy E above V(r e d ge ), 
in which case the orbits will generally be hyperbolic-like while keeping the self-similiarity of the configuration. 
A radial orbit with a large apocenter at r apo would have a period 
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y = ^- (30) 

Tedge 



where we have neglected the finite size b and Newtonian potential, and rescaled all length with the radius r apo and 
Tedge- Note that an radial orbit which just escapes with E = V(r e d ge ) will take infinite amount of time to reach r ec [ ge 
because of the linear decline of the radial speed dr/dt oc (r e dge — f) near the edge. 

Escape speed, the edge of potential, and the cutoff of the effective dark halo 

Following [23], we can define the effective DM mass Medm = (g — 9n)t 2 /G. This is the equivalent Newtonian mass 
of the DM to generate the same potential as in MOND. Neglecting the finite size b, we find Medm has a positive 
part linear to r and a negative r 3 part, hence the corresponding effective DM density Pedm = d ^r°dr nas a positive 
1/r 2 part , and a uniform negative part p c — ?>P c /c 2 = —2p c with the point of zero effective density at 

r cut = ^ = (GM«o) 1/4 ^. (31) 

At this radius the circular speed v C i r (r cut ) and the total dynamical mass inside can be estimated by M cut — 
v1i r T C utG~ l — M + m + 7^\/GMaor cu t. A circular orbit bound at r = r cut has a period T cir — \f2ixT. For ex- 
ample the cut off radius of the Milky Way is about 2000 kpc, where the peak of the effective potential curve is (see 
Figl). 

A working definition of an escaping orbit is an orbit which reaches r cut , since a radial orbit which reaches an 
apocenter r apo — r cut would have a period T ra( n a i ~ (1-5 — 1.7)r ~ 16 — 19Gyrs. Such an orbit takes slightly too long 
to return over the Hubble time 14Gyrs. So the escape speed at any radius r is defined as 

v esc (r) = ^2V(r cut ) - 2V(r)\ j=0 . (32) 

E.g., we can estimate the escape speed near the solar neighborhood r = 8kpc adopting r cut ~ 1500kpc and 
(GMcto) 1 / 4 = 180 km/s for the Milky Way. Our analytical formulae predicts a local escape speed v esc ~ 



1/4 



1/4 



y/2 In Si x (GMa ) 1/4: km/s - 580 (%f\ -580 % km/s, consistent with the observed value ~ 550 km/s [23] 
for solar neighbourhood stars. The so-called external field effect is not as critical here as in [23], but our prediction 
depends on the cosmological constant and depends on the satellite-galaxy mass ratio m/M through the function || , 
which is unity for a single star orbiting a galaxy; interestingly the escape speed is smaller for a more massive satellite 
in MOND, however, the difference is fairly small as long as the satellite is less massive than the mass of the host 
galaxy, (see Fig. 1). 
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FIG. 2: shows the orbits computed for the LMC-MW pair with a LMC baryonic mass M/5 (small circles) or M/10 (lines) 
respectively for the past 12 Gyrs; the LMC is presently moving to the left with 380km/s at 45kpc below the MW, which is 
held at the origin (shown by the red dot). A generic difference with Newtonian Keplerian orbits is that the MONDian orbits 
depends on the mass ratio, and the direction of the pericenter or apocenter changes as the nearly elliptical orbit precess. 

TIMING THE MAGELLANIC CLOUD, THE ANDROMEDA, THE BULLET CLUSTER 



Timing is a classical argument which assumes two presently neighbouring bodies were very close at birth within an 
expanding Hubble flow away from each other. The mutual gravity eventually overcomes the Hubble flow and brings 
their orbits close to each other again. Their orbital or dynamical age must then be close to the Hubble time. A few 
timing applications of our analytical result is shown in Fig. 2, Fig. 3 and Fig. 4. The relative distances of these are 
such that the LMC moves inside the Milky Way, and is presently near the pericentre, and the M31 is approaching us 
from about 800 kpc after falling out of the Hubble expansion. 

First consider the dynamical age of the Large Magellanic Cloud (LMC), a satellite on an non-circular orbit around 
the galaxy. Recent observations found the LMC moves with a speed of 360 km/s almost tangentially at a distance 
45-50 kpc from the Milky Way, which is too fast for a Newtonian Dark Halo to keep it bound[23]. In our scenario the 
MONDian Milky Way has a baryonic mass M\ = 5 x 10 10 solar masses, and the LMC has a baryonic mass M^ = Mi/5, 
hence (GaoM) 1 / 4 = 180km/s, we find r cut = 1500 kpc, r max = exp((360/180) 2 /2) x 45kpc = 350kpc = r cut /5. We 
estimate the oscillation period in the radial direction T rac iiai = 0.3r = 4 Gyrs. The period in the circular direction 
is about 1.4 times longer. This implies that the Large Magellanic Cloud could be bound and have enough time to 
circulate the Milky Way more than once in MOND (cf. Fig. 2). This is consistent with the standard scenario where 
the Magellanic stream is pulled out from the Magellanic Clouds during one of its pericentric passages of the Milky 
Way, and possibly interactions among themselves [11]. We also note that the baryonic mass of the Milky Way 
adopted for modeling the LMC is roughly consistent with the local escape speed and the orbit of the Sgr dwarf. As 
shown in Fig. 1, in particular our model orbit has a pericenter and apocenter of 10 and 50 kpc respectively, and a 
present 280 km/s tangential velocity at a radius 16 kpc from the Milky Way. These parameters are consistent with 
the observational data and typical orbits found in previous dark halo models [13]. 

We have also computed the orbit of the M31, whose parameters are not well-determined observationally. The 
M31 galaxy is presently at ~ 800kpc and is moving at 130km/s towards us. Its tangential velocity is unknown, 
although perhaps small. We study two models: in one model M31 has an equal baryonic mass as the Milky Way 
M<x = M\ =5x 10 10 solar masses and M31 has a nominal small transverse velocity of 40km/s. In another model 
M31 is twice as heavy as the Milky Way[5] and has a larger transverse velocity lOOkm/s to avoid collision with the 
M33[12]. 

Shown in Fig. 3 are the orbits in the past 12 Gyrs, earlier than which the disks of the Milky Way and M31 are 
likely not yet formed. MOND prefers either a small baryonic mass for the Local Group, or a significant transverse 
velocity of the M31, something that can be falsified by future measurements of the transverse velocity. A transverse 
velocity helps to keep the pericenter distance large, about 200-300kpc, hence any tidal effect from the M31 on the 
LMC is small, which is perhaps a desirable feature to bound the LMC within 45-450 kpc of the Milky Way. However, 
if we adopted a radial orbit and a large mass for the M31 as in the standard interpretation of Local Group timing 
with a Keplerian orbit of M31 and MW dark halos [4, 18], there would be some tension between the long age of the 
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FIG. 3: Shown are two possible orbits for M31 (green)-MW (red) for the past 12 Gyr: the left panel shows a possible nearly 
radial orbit with the M31's baryonic mass the same as that of the Milky Way (M = 5 x 10 10 solar masses) and the right panel 
shows a significantly non-radial orbit of the M31 with a mass 2M. The origin is set at the present position of the MW, while 
the M31 is presently 800 kpc away on the right. The arrows indicate the force towards the center of mass. 
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FIG. 4: shows the orbit of the Bullet subcluster with respect to the main cluster (green dot fixed at origin) with a mass ratio 
1:2 (red circles) and 1:3 (red line) respectively for 9 Gyrs earlier than its present age; the Bullet is at 700kpc to the right of 
the main cluster and is moving to the right at 4000km/s; both models assume the same combined mass Mi + M2 = 6 x 10 14 
solar masses. 



universe and short orbital period of M31-MW binary predicted in MOND, implying that the two systems had an 
earlier flyby, and are coming close to each other for the second time. Our model cannot yet include any acceleration 
of the whole M31-MW binary towards the Virgo cluster, which creates the so-called external field effect, which can 
generally reduce the MOND force or potential, and lengthen the period [8, 23]. 

To time the Bullet Cluster, which is at z = 0.3 when the age of the universe is 10 Gyrs, we can set T rat n a i ~ 10 
Gyrs. A possible orbit is shown in Fig. 4. Allowing for some hydrodynamical effect [21], we set the speed of 

encounter J 2 In ^(GMa ) 1/4 - (3000 - 4000) km/s at the present separation r = 700 kpc, we find (GMo ) 1/4 > 

3000/-y/21n(7) ~ 1500 km/s, hence the matter M > 2.5 x 10 14 solar masses, much larger than the combined baryonic 
content < 10 14 solar mass for both systems. This implies the need for non-baryonic matter, perhaps (sterile) neutrinos 
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in both systems [19]. These crude estimates are consistent with the previous findings [1, 17]. 

SUMMARY 

In short, previous tests of MOND have largely been limited to fitting rotation curves or velocity dispersion curves 
inside axisymmetric galaxies or galaxy clusters. A few studies in the literature on non-axisymmetric configurations 
have relied primarily on numerical codes [6, 9, 16, 22, 23]. The numerical complexity severely limits our theoretical 
intuitions on this class of non-linear theory of gravity. 

Here we derive the modified Kepler's law analytically for a two-body system and for restricted many-body problem 
in the context of two versions of theories for MOND. We demonstrate the powerful use of modified Kepler's law by 
applying our analytical results to make predictions of the orbital motions for real systems. These analytical results 
are also useful for testing numerical codes (e.g., [6]) and for getting intuitions. It appears that in the case of the 
Bullet Cluster a pure MOND without non-baryonic matter (e.g., neutrinos) is not enough to explain the fast motions 
of the bullet. On the other hand the timing of the M31's orbit would imply an uncomfortably low mass for the M31 
unless its orbit towards the Milky Way has significant amount of angular momentum, which is a strong prediction 
for MOND to survive in the context of the two-body problem in the Local Group. The baryonic mass of the Milky 
Way ~ 5 x 10 10 solar masses seems enough to consistently explain the rotation curve [24], the local escape speed, and 
the morphology and kinematics of the Sgr stream (Fig. 1); note in our approximation we have neglected the effects of 
the detailed form the MOND /i function in these regimes. The orbit of the LMC is bound around the Milky Way in 
MOND, however, with only two pericentric passages in the past 12 Gyrs, so it remains to be seen if these pericentric 
passages are enough to generate the detailed morphology of the Magellanic Stream. 

To conclude, our analytical formulae for the MOND gravity in two-bodies provide some tools for studying MOND 
beyond rotation curve fitting, where MOND has been very successful. We propose to use the timing argument in two- 
body problem to probe the dynamics of the LMC, M31 and the Bullet clusters in MOND. The ultimate falsification 
or proof rests on more detailed numerical modeling with improved kinematic data. 

Appendix A: Motion-independent force and finite-size correction 

Assume the body mi is a Kuzmin-Hcrnquist disk-bulge flattened system (as introduced in Shan et al. 2008) [20] 
with two imaginary centers at | ± k| above or below the plane of the disk [4], one can apply the formulae as if mi is 
a spherical Hcrnquist body of scale length b centered on a point below the plane whenever the body m 2 is above the 
plane, and vice versa. So the equation of motion for mi and m 2 are 



m 2 d 2 r 2 mid 2 ri d Gm\m 2 2y / Gao(mi + m 2 ) 3 

dt 2 alt 2 =+ 9r^|r 2 -ri±k|+6 3 

For two spherical particles with k — 0, the mutual force 




3/2N 



Gm 1 m 2 -v /G ( TO i + TO 2) 3 »o - 2 / ^ / m l » • » 

F 12 = mioi = m 2 a 2 = —r^ + ~ " T" T » - = o I 1 ~ >. T > ( 34 

(ri2+by n 2 + b Sy ^—^\mi+m 2 J J 

where we have opted to smooth our point-like particle with a common scale b by a Hcrnquist-smoothing Kernel, which 
is not always rigorous, but allows a non-divergent estimation of forces in situations where the bodies are overlapping. 
The essential thing is to keep the forces rigourous for point mass and that F\ 2 — F 2 \ in general such that the system 
conserves total momentum. The forces on the two bodies are in opposite directions, keeping their center of mass fixed. 

Appendix B: Alternative expression for M in binary configuration, and in symmetric many-body 

configurations 

Following expressions for the effective MONDian mass M have better asymptotic behavior for N — 1 binary 
configuration with two masses M and m, and for a central mass M plus a N > 1 symmetric identical particles of mass 



m/N. 
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0)l / fl/2 rn l/2 



M 3 / 2 



2M(M + m) + (2/3)m 2 2m 1 / 2 
(M + to) 3 / 2 + M 3 / 2 3N 1 / 2 



,3/2 



2(M + to) 3 / 2 



otherwise 



if AT = 1 



(35) 



For equal mass binary m = M, and TV = 1, we have M = 2 [l — |] [| + 575] to ~ 1.2m. 

Appendix C: expressions for mjv for N > 1 identical masses distributed on a regular polyhedra 



It is somewhat lengthy but straightforward to calculate the Newtonian force vectors between N particles distributed 
on regular polygons and polyhedra of radius r. Summing up the forces on each particle, one can find the forces 
are indeed centripedal, the acceleration is given by GrfiN/r, The expressions for a^ = ^^- are found as follows: 



v'2 



for a N — 6 points on octahedron and ag, 



a 6 24 ' 6 _ 

totally 5 possible regular polyhedra or platonic solids, with N 
cube, icosahedron, and dodecahedron. 



l+ A /372+v / 3 



.32 



for a N 



points on a cube. There are 



4,6,8,12,20 vertices for tetrahedron, octahedron, 
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